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We extend the symmetrized density matrix renormaliza- 
tion group (SDMRG) method to compute the dynamic non- 
linear optic coefficients for long chains. By computing correc- 
tion vectors in the appropriate symmetry subspace we obtain 
the dynamic polarizabilities, ay(w), and third-order polar- 
izabilities jijki (w, ui, w) of the Hubbard and "U — V" chains 
in an all — trans polyacetylene geometry with and without 
dimerization. We rationalize the behavior of a and j on the 
basis of the low-lying excitation gaps in the system. This is 
the first study of the dynamics of a fermionic system within 
the DMRG framework. 

PACS Numbers: 42.65.Sf, 61.82. Pv and 78.20.Bh 



The dynamic linear and nonlinear response functions 
of finite interacting systems have most commonly been 
obtained from an explicit computation of the eigenstates 
of the Hamiltonian and the matrix elements of the ap- 
propriate operators in the basis of these eigenstates |jj . 
This has been the most widely used method, particu- 
larly in the computation of the dynamic Nonliear Optic 
(NLO) coefficients of molecular systems and is known as 
the sum-over-states (SOS) method. In the case of model 
Hamiltonians, the method that has been widely used to 
study dynamics is the Lanczos method The spec- 

tral intensity corresponding to an operator O is given by 
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where \G > is the ground state eigenvector of the Hamil- 
tonian with eigenvalue Eq, z = lj + ie, to is the frequency 
at which the response is sought, e is the mean life time 
parameter and H is the Hamiltonian. In the Lanczos 
method, I(uj) is computed as a continued fraction, 
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wherein the coefficients a$ and hi are respectively the 
diagonal and the off-diagonal matrix elements in the 
tridiagonal matrix representation of the Hamiltonian ob- 
tained in the Lanczos procedure. In this method, there 



is an implicit truncation of the Hilbert space due to the 
smaller size of the tridiagonal matrix compared with the 
total dimensionality of the space spanned by the Hamilto- 
nian [gj. Therefore, the dynamic quantities computed by 
this technique are approximate even though the ground 
state obtained is exact. 

In the context of NLO properties of the Hubbard 
and extended Hubbard Models it was shown by Soos and 
Ramasesha that the model exact dynamical NLO coeffi- 
cents could be obtained by solving for correction vectors 
. If we define the correction vector <fP^ (uj) by the equa- 
tion 

{H-E n -z)^ 1 \u J ) = -d\G>, (3) 
then the spectral function, I(uj), can be expressed as 

J(w) = —Im < G\6 ] \<t> {1) {u) > (4) 



The correction vector is solved for in the basis of 
the configuartion functions, which is also the basis in 
which the Hamiltonian matrix is set-up for obtaining the 
ground state. Given the ground state and the correc- 
tion vector, it is straightforward to compute the spectral 
function. This method is quite general and has been em- 
ployed in the computations of dynamic NLO coefficients 
of a wide variety of Hamiltonians || . The inhomogeneous 
linear algebraic equations encountered in this method in- 
volve large sparse matrices and an iterative small matrix 
algorithm, which runs parallel to the Davidson algorithm 
for eigenvalue problems, gives rapid convergence for the 
solution of the system of equations |J . 

The two correction vectors 4>f^ (bj\) and <fy? (a;x,W2) 
encountered in the computation of polarizability and 
third-order polarizability are solved for from the follow- 
ing linear equations . 
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where fliS are the dipole displacement matrices and other 
quantities are as defined in eqn.(l). The DMRG method 
PIPI as implemented, readily provides us with the 
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ground state and the Hamiltonian matrix. The matri- 
ces of the dipole operators are constructed in the DMRG 
scheme by renormalizing the matrix representations of 
the dipole operator corresponding to the left and right 
parts of the system using the density matrix eigenvector 
basis in a way which is completely analogous to the corre- 
sponding Hamiltonian operators for the fragments. The 
matrix representation of the dipole operators for the full 
system are obtained as direct products of the fragment 
matrices analogous to the way by which the full Hamil- 
tonian matrix is constructed. The dipole displacement 
matrices are obtained by subtracting the corresponding 
components of the ground state dipole moments from the 
diagonal elements of the dipole matrices. 

In terms of these correction vectors, the components of 
the polarizabilities, a^, and third-order polarizabilitics, 
Jijki , can be written as 

Oij(u) =< ^MI^-IG > + < ^(-wJIftlG >, (7) 



. (8) 

where the operator P generates all the permutations: 
(— uj a , i), (u>x,j), (u>2, k) and (a>3, 1) leading to 24 terms for 
lijki with w a = —lo\ — u>2 — W3. The tumbling averaged 
a and 7 can be defined as 

1 3 1 3 

«=3^ a " ; ^ = YsX! ( 2 >« + 'ftui) ( 9 ) 
i=l ij = i 

to allow comparison of the calculated NLO response with 
experiments on systems containing molecules in random 
orientations 

The DMRG method as usually implemented, does not 
exploit all the symmetries of the system. In the case 
of model Hamiltonians for polymers, the system poss- 
eses total spin symmetry, reflection symmetry about the 
middle of the chain and in some cases the electron-hole 
symmetry. These symmetries ensure that a given cor- 
rection vector spans only a symmetrized subspace and 
not the entire Hilbert space of the Hamiltonian of the 
given system. The correction vector lies in the symmetry 
subspace which is connected to the ground state by the 
electric dipole operator. For lo values corresponding to 
resonance between the ground state and eigenstates in 
that particular symmetry subspace, Ihs of eqns(5) and 
(6) become singular for e = and present numerical dif- 
ficulties even when solving them for reasonable nonzero 
e values. However, if we do not exploit the symmetries of 
the Hamiltonian, we encounter sigularities in eqns(5) and 
(6) even for those ui values corresponding to eigenstates 
of the Hamiltonian found in other symmetry subspaces 
which are not connected to the ground state by the dipole 



operator. Therefore, numerically it would be impossible 
to obtain the correction vectors at these frequencies and 
thereby the associated response of the system. 

For example, in Hubbard chains at intermediate cor- 
relation strengths, a triplet excited state lies below the 
lowest singlet state in the ionic B subspace ]l3| |. The 
states in the ionic B space are connected to the ground 
state via one-photon transitions. The resonances in po- 
larizability are thus expected only at frequencies corre- 
sponding to the energy levels in the ionic B space, relative 
to the ground state. However, we can not solve for the 
correction vector using equation (5) at an excitation en- 
ergy corresponding to the the energy of the lowest triplet 
state. Thus the technique of correction vectors will not 
be able to give the complete dispersion of the polariz- 
abilities upto the first one-photon resonance, unless in- 
terferences due to spurious intruders such as the triplet 
states are eliminated by suitably block-diagonalizing the 
Hamiltonian matrix. This problem of intruders becomes 
more severe with increasing system size due to increas- 
ing number of intruder states lying below the frequency 
corresponding to the first "true" resonance. 

We have exploited the electron-hole symmetry, the re- 
flection symmetry of the polymers and spin parity to 
block-diagonalise the Hamiltonian p"4j . The electron- 
hole symmetry divides the Hilbert space into ionic and 
covalcnt spaces. The ground state is in the covalent space 
while all the dipole allowed excitations from the ground 
state lie in the ionic space. Use of parity conservation 
divides the Hilbert space of the Hamiltonian into even 
and odd parity spaces corresponding to even and odd to- 
tal spin states. We have considered a single excitation 
frequency in all our calculations (lo = u>i — u>2 = L03) and 
have exploited sparseness of all the matrices to improve 
upon the computational efficiency. We have employed 
the finite size DMRG algorithm in some cases to check 
the convergence of the DMRG results. In these cases, 
the spatial symmetry is exploited only at the end of the 
DMRG procedure when the left and the right density 
matrices correspond to fragments of the same size, i.e., 
at the end of each finite size iteration. 

We have computed the dynamic linear polarizabilitics 
(a(u))) and third-order polarizabilities (7(^,^,0;)) cor- 
responding to the third harmonic generation (THG), for 
the Hubbard and "U—V" chains of upto 20 sites with and 
without dimerizations, for many values of the parameters 
of the Hamiltonian, 

H= [l-^(-l) i ][-*(aJ )CT aj,a+aj i(T ai, CT )]+C7n i0 .ri i _ .+ 

<ij>,cr 
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where a\ a (ai^),hi a , t and U have their usual meaning. 
V is the nearest neighbour interaction parameter and the 
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term in V is nonzero only when the nearest neighbours in 
question are charged, in which case, the interaction is re- 
pulsive for like charges and attractive for unlike charges. 
8 is the dimerization parameter of the system. The ge- 
ometry of the Hubbard and " U — V" chains required for 
the computation of the NLO coefficients is chosen to cor- 
respond to all — trans polyacetylene configuration, with 
a bond angle of 120° and bondlength of 1A, for the uni- 
form chain. Dimerization (8 > 0) leads to proportionate 
alternation in bond lengths as well as alternation in the 
t and V parameters. All the computations have been 
carried out at a single frequency corresponding to an ex- 
citation energy of O.leV^ and the life-time e is chosen to 
be 0.00W. 

We have compared the DMRG results with a cut-off 
of to = 200 (retaining the dominant 200 density matrix 
eigenvectors in the DMRG scheme) with the model ex- 
act a and 7 values obtained from the correction 
vector method for chains of upto 12 sites. For the 8-site 
problem, the DMRG calculation is exact for this cut- 
off and the DMRG results compare with exact results 
to numerical accuracy. In the case of the 12 site uni- 
form Hubbard chain with U/t = 4, the model exact a is 
5.343 x 10" 24 esu while the DMRG a is 5.293 x 10~ 24 
esu. The model exact 7 value is 598.3 x 10~ 36 esu while 
the DMRG 7 value is 589.3 x 10 -36 esu. The error in the 
dominant a and 7 components namely, a xx and "f X xxx is 
much smaller. The model exact a xx is 14.83 x 10~ 24 esu 
while the DMRG a xx is 14.79 x 10~ 24 esu. As for 7, the 
model exact "/ xxxx is 28 73 x 10 -36 esu while the DMRG 
j xxxx is 2870 x 10~ 36 esu. 

In figs.l, we dispay the dependence of a on the chain 
length for different values of U/t, for uniform (fig. la) 
and dimerized (fig. lb) Hubbard chains. The polarizabil- 
ity decreases with increasing correlation strength in both 
the dimerized and the uniform chains. For the same chain 
length the uniform chains have, as expected higher 
polarizability than the corresponding dimerized chain at 
every value of U /t we have studied. For the uniform 
chain, the average polarizability, a, for weak correlations 
exhibits a nice power law dependence on chain length 
with an exponent of 2.022 ± 0.001. However, for stronger 
correlations, the polarizability deviates from a power law 
dependence and seems to show a size-consistent varia- 
tion at longer chain lengths. The chain length at which 
the the change over from the power law behavior occurs, 
systemetically reduces with increasing U ft. As the chain 
dimerizes, the range over which power law behavior is 
observed decreases (Fig. lb). Eventually, in the limit 
8 = 1.0, we would observe only size-consistent depen- 
dence as, in this limit the chain breaks down to nonin- 
teracting dimcrs. 

In figs. 2, we show the log-log plot of average third- 
order polarizabilty 7 with the chain length for chains of 
upto 20 sites with (fig. 2a) zero and (fig. 2b) nonzero 
bond alternation. The variation of 7 with chainlength 



is similar to that found for a both in the dimerized and 
undimerized cases. Nonzero 8 leads to a decrease in the 
7 value and the decrease is smaller for higher U/t values. 
An exponent of 5.570 ± 0.001 is found for 8 = 0.0 and 
U/t = 2.0, while the exponent is 4.00 ± 0.01 for 8 = 0.09 
and U/t = 2.0. In both these cases, the power-law behav- 
ior is observed upto the maximum chain length we have 
studied. It is interesting to note the strong dependence 
of the exponent on the dimerization parameter. 

In the " U — V" model the dependence of the polariz- 
ability on the chain length is very simialr to that in the 
Hubbard model. However, the "U — V" model is more 
polarizable, for the chosen model parameters and fre- 
quency. While the dimerization, 8, decreases the polar- 
izability of the chain, the nearest neighbour interaction, 
V, increases the same. The exponents for a, where the 
power law holds (U/t = 2.0, V/t = 1.0) is 2.320 ± 0.001 
for the uniform chain and 1.64 ± 0.01 for the dimerized 
chain(<5 = 0.09). 

In figs. 3 we present the dependence of 7 on chain 
length for various values of U/t for the uniform (fig. 3a) 
and dimerized (fig. 3b) chains, in the "U — V" model 
with V/t fixed at 1.0. There are some very interesting 
differences in the variation of 7 between the Hubbard 
and the "U — V" chains. The dependence of the 7 of 
the two chains on U are similar for U > 2V. In this 
the SDW regime, the third-order polarizability in the 
" U — V" model is larger than that in the Hubbard model 
for corresponding chain lengths both with and without 
dimerization. However, when U = 2V, which is the 
crossover point from the SDW to the CDW regime, we 
find that the Hubbard chains have higher third-order po- 
larizability than their counterparts in the " U — V" model. 
This trend is observed only for 7 and is not seen for the 
polarizability a. To understand this behavior we have 
studied the variation of the lowest one — photon gap (fig. 
4a) and the lowest two — photon gap (fig. 4b) as a func- 
tion of the chain length in the Hubbard model and the 
"U-V" model at the SDW/CDW transition point. The 
gap to the lowest one-photon state in the former is higher 
than that in the latter at all chain lengths, independent of 
8. Therefore the polarizabilities of the " U — V" model at 
the transition point are higher than that of the Hubbard 
model. The third-order polarizability also has contribu- 
tions from the two-photon states. The lowest two-photon 
state in the " U — V" model for these parameters is at a 
higher energy than in the Hubbard model for all the chain 
lengths and 8. Consequently, the two-photon contribu- 
tion to 7 is higher for the Hubbard model than for the 
"U-V" model. 

For values of V/t > U/2t (CDW regime), the opti- 
cal gap decreases sharply and the chosen excitation fre- 
quency, u>, of O.ley is above the first three-photon res- 
onance in the system. In this regime, 7 is also found to 
have a negative sign for sufficiently long chains (N > 10 
sites). 
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To conclude, we have demonstrated how dynamic NLO 
responses can be obtained within the DMRG procedure. 
We have applied the method to Hubbard and " U — V" 
chains in the CDW and SDW regimes. 
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Fig.l 

Log-Log plot of the average polarizability (a) in 10~ 24 
esu versus chain length L in A for Hubbard chains of 
upto 20 sites in all — trans polyacetylene geometry for 
(a) 5 = and (b) S = 0.09, for three different values of 
U/t. 
Fig.2 

Plot of the log of average third-order polarizability (7) 
in 10~ 36 esu versus log of the chain length L in A, for 
Hubbard chains, with three different values of U/t for (a) 
5 = and (b) 5 = 0.09. 
Fig.3 

Log-Log plot of the average (7) in 10 -36 esu versus chain 
length L in A for "U — V" chains for three values of U/t 
and for V/t = 1 for (a) 5 = and (b) 5 = 0.09. 
Fig.4 

Dependence of (a) one-photon gap (ground state to 
l 1 !?") and (b) two-photon gap ( ground state to 2 1 j4+) 
on l/N. N is the number of sites in the chain. 
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